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Abstract. We present a functional renormalization group analysis of superconduc- 
tivity in the ground state of the attractive Hubbard model on a square lattice. Spon- 
taneous symmetry breaking is treated in a purely fermionic setting via anomalous 
propagators and anomalous effective interactions. In addition to the anomalous in- 
teractions arising already in the reduced BCS model, effective interactions with three 
incoming legs and one outgoing leg (and vice versa) occur. We accomplish their inte- 
gration into the usual diagrammatic formalism by introducing a Nambu matrix for the 
effective interactions. From a random-phase approximation generalized through use of 
this matrix we conclude that the impact of the 3-1-1 effective interactions is limited, 
especially considering the effective interactions important for the determination of the 
order parameter. The exact hierarchy of flow equations for one-particle irreducible 
vertex functions is truncated on the two-particle level, with higher-order self-energy 
corrections included in a scheme proposed by Katanin. Using a parametrization of 
effective interactions by patches in momentum space, the flow equations can be inte- 
grated numerically to the lowest scales without encountering divergences. Momentum- 
shell as well as interaction- flow cutoff functions are used, including a small external 
field or a large external field and a counterterm, respectively. Both approaches pro- 
duce momentum-resolved order parameter values directly from the microscopic model. 
The size of the superconducting gap is in reasonable agreement with expectations from 
other studies. 
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1. Introduction 

Renormalization group (RG) techniques facilitate the systematic investigation of 
interacting many-particle systems even in situations where perturbation theory breaks 
down. A central idea underlying many renormalization group methods, pioneered by 
Wilson m, is to take into account the degrees of freedom successively, ordered by 
their kinetic energy. Formally, this is accomplished by integrating a set of differential 
equations called flow equations, where some energy scale (usually a cutoff) plays the 
role of the flow parameter. The RG flow generates a family of effective actions which 
interpolates smoothly between the bare action of the system and the final effective action 
where all degrees of freedom have been integrated out. 

The first significant applications of RG methods to interacting Fermi systems in 
more than one dimension were undertaken by mathematicians aiming at a rigorous 
control of certain systems and their properties [2|, [3|, H] . The usefulness of RG concepts 
for Fermi systems was soon recognized also by physicists [HIE]. RG calculations for Fermi 
systems are often termed "functional" as they involve the flow of functions instead of 
just a few running couplings. An early computational success of the functional RG 
(fRG) was unambiguous evidence for d-wave superconductivity in the repulsive two- 
dimensional Hubbard model [H [3 [9] . 

Of particular interest in interacting Fermi systems is the spontaneous breaking of 
symmetries. Formally, such a symmetry breaking occurs if a product of two fermionic 
fields - a bilinear - which does not appear in the original Hamiltonian acquires a nonzero 
expectation value. This bilinear can be replaced by a bosonic field through a Hubbard- 
Stratonovich transformation, which simplifies the treatment of the order parameter and 
its fluctuations (lOj. However, these simplifications come at the price of introducing a 
bias towards breaking the symmetry in the channel fixed by the Hubbard-Stratonovich 
field. Working with the original fermionic degrees of freedom this bias can be avoided. 
However, the symmetry-broken phase is not easily accessible by a purely fermionic RG 
flow. Calculations starting from the symmetric phase run into divergences at a finite 
energy scale, that is before all degrees of freedom have been integrated out. 

Most previous functional RG studies of interacting Fermi systems are based on a 
flow equation for the generator of the one-particle irreducible (IPI) vertex functions 
This equation is equivalent to an infinite hierarchy of flow equations for the IPI 
vertex functions. Truncating the hierarchy in a weak coupling expansion, the RG flow 
has been computed for various low- dimensional Fermi systems [8], [121 [131 [131 [13 [16]. 
Instead of the common choice of a momentum cutoff, it was found to be advantageous in 
certain situations to use as flow parameter the physical temperature [18], the strength 
of the interaction [19], or a cutoff in frequency space [201 [2T] . 

Spontaneous symmetry breaking within the functional RG framework has 
frequently been studied by introducing bosonic fields for the order parameter, leading to 
a coupled flow of bosons and fermions [22l [23l [24]. A simple way to access spontaneous 
symmetry breaking in a purely fermionic setting is to combine the fRG with a mean- 
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field calculation for the low-energy degrees of freedom [TTj. However, order parameter 
fluctuations are thereby discarded. To continue the fermionic functional RG flow into 
the symmetry-broken phase, an improved truncation of the hierarchy with an efficient 
treatment of self-energy corrections, first formulated by Katanin [25] , is very promising. 
Employing this improvement, the Bardeen-Cooper-Schrieffer (BCS) mean-field model 
can be solved exactly already by a low-order truncation [26], if a small external field 
is used to trigger the symmetry-breaking. The vaunted divergence is then regularized, 
but there is still a flow regime with strong coupling. The method was extended to the 
breaking of a discrete symmetry and finite temperatures [27] . Employing the interaction 
flow [19], the external field can be compensated by a counterterm, thus enabling the 
study of first-order phase transitions which were previously out of the method's reach 
[28] . Furthermore, this approach yields results at zero external field without having 
to pass through a strong coupling regime in the case of discrete-symmetry breaking. 
Strong attractors for the flows of various counterterm strengths were found. 

In this work, we apply the fermionic fRG to the attractive Hubbard model, using 
both the external-field as well as the counterterm approach. The repulsive Hubbard 
model [291 [30l [31] was originally introduced to study ferromagnetism of itinerant 
electrons, but has since been employed to elucidate many aspects of strong correlations 
in solid-state physics. The attractive version exhibits s-wave superconductivity under 
certain conditions, including the case we study here, namely for a square lattice at 
zero temperature and weak or moderate coupling [32]. It is a good testbed because 
the correct order parameter is affected by fluctuations even in the weak coupling limit 
[33] , while the model is fairly well-studied, allowing comparisons of the results obtained 
using the new method. Lattice fermions with attractive interactions can nowadays be 
reahzed by cold atoms in optical lattices, and s-wave superfluidity in such systems has 
been observed very recently [3^ . 

Our goals in this work are threefold. First, we would like to study the role of 
anomalous effective interactions with three ingoing particles and one outgoing particle 
(or vice versa), which arise in systems with both Cooper and forward scattering. We 
find that these interactions have a limited impact. Second, we would like to determine 
whether the application of the Katanin-truncated fRG to the Hubbard model works as 
well as for mean-field models. For the momentum-shell approach with a small external 
field, we obtain reasonable flows resulting in order parameter values comparable to 
the literature. It is possible to choose the strength of the external field two orders 
of magnitude below the final order parameter value without running into artificial 
divergences. For the interaction-flow approach, we find remnants of the strong-attractor 
behavior known from the mean-field models. Third, we compute the superconducting 
order parameter in the ground state of the attractive Hubbard model as a function of 
coupling and filling. 

This paper is structured as follows. In section [21 the attractive Hubbard model 
is introduced formally. The formalism permitting the study of symmetry breaking 
in general and the 3 + 1 anomalous effective interactions in particular is detailed in 
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section [31 It is employed in the context of a generalized random-phase approximation 
(RPA) of the model's behavior in section |H Section [5] presents the result of the fRG 
studies, containing momentum-shell and interaction flows as well as a study of the order 
parameter strength and a comparison to earlier calculations. 

2. Attractive Hubbard model 

We discuss the attractive Hubbard model on a two-dimensional square lattice specified 
by the Hamiltonian 

ks kik2k:i 

N denotes the number of lattice points, ?7o > is the onsite attraction, and k takes all 
values in the first Brillouin zone of the square lattice. The dispersion relation 

^fc = — 2t(cos kx + cos ky) — At' cos k^ cos ky — fi (2) 

contains a probability amplitude t for nearest-neighbor hopping, a probability amplitude 
t' for next-nearest-neighbor hopping, and a chemical potential /i. An external field 

J2 {\.t{k)cl^icU + KAk)c^^c_^^) (3) 

k 

is added to H to trigger the instability towards superconductivity. This field is either 
kept small or turned off gradually in the fRG flow to approximate the situation with 
spontaneous symmetry breaking. 

3. Nambu formalism 

3.1. Functional integral and field substitution 

We study the partition function 2 in a functional integral representation 

Z = J V{^,^)e^^^''f'\ (4) 
within the Matsubara formalism, employing the grand canonical action 

S = ^{iuJn - ^k)lpks'^ks - ^ (Acxt(fc)^-fci^feT + 
ks k 

+ Uo ^fclT^fc3T^fc2i^fcl+fc2-fc3i (5) 

ki,k2M 

where ip and are the Grassmann fields matching c and c\ respectively. Furthermore, 
k = [ioOn, k) and volume as well as temperature factors are implicitly contained in the 
summation symbol a convention we also adopt in the following. We employ the 
substitution 



i^kh (pk+ = i'kh <Pk- = ^-ki, 4>k- = i'-ki (6) 
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to obtain 




kT_-k2-4>k:i + 4>-k2- 




(7) 



kik2k3 



3.2. Propagator and effective interaction 



In (I7j), the external field appears on the off-diagonal of a matrix specifying the action's 
quadratic part. According to the Dyson equation, inverting this matrix after subtracting 
the self-energy leads to the propagator or one-particle Green's function 



-lUOr 



ujl + E{kY 



A*{k) 





(8) 



(9) 



where E{k) = ^ {^k — + I^P and is the normal self-energy. A denotes the 

superconducting order parameter and contains the anomalous self-energy as well as the 
external field. To simplify the diagrammatics, we antisymmetrize the interaction part 

Uo 5fci+fc2,fc3+fc4 {^Ni + Sn2~ ~ ^N2 + ^Ni~) {^Na+^Ni- — Sn4+^N3-) 01020304, (10) 



1...4 



where i, if used as an index, is short for kiNi and 2 and 4 are exchanged in comparison 
with ([7j). If (NiN^) labels the rows and {N2N4) the columns of a matrix to the Nambu 



base {(- 



-)}, the antisymmetrized bare interaction reads 



Uo = f/o 



(0 








1 \ 








-1 








-1 
















J 



4i 



+k2,ks+k4- 



(11) 



By employing the Nambu field substitution ([6]), a formalism containing four potentially 
different normal (particle-number conserving) and twelve potentially different anomalous 
(particle-number non-conserving) effective interactions arises. Here, as we concentrate 
on conventional s-wave pairing, we require time reversal invariance to hold. This implies 
the following behavior of an effective interaction component under flipping of all its 
Nambu indices: if an odd number of the Nambu indices are (+), the sign of the effective 
interaction's real part is flipped; otherwise, the sign of the imaginary part is flipped. 
Exploiting this symmetry, we assign symbols to the various effective interactions: 



U 



X* 

u* 



V* 



V 

w 



u \ 

ill 



(12) 
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The f2j are collectively referred to as 3 + 1 effective interactions as they have three 
incoming and one outgoing leg or vice versa if drawn according to the ip-field 
representation of ([5]). They are generated in the symmetry-broken phase by diagrams 
of the type in figure [H is a 4 + anomalous effective interaction arising only in the 




Figure 1. Diagram in ^/i- field language generating an effective interaction with three 
incoming arrows in the symmetry-broken phase. Bold lines with one (two) arrow(s) 
represent full normal (anomalous) propagators, wiggiy lines represent bare interactions. 

symmetry-broken phase. V, X, and U are normal 2 + 2 effective interactions. Only V 
and U have non-zero bare values. 

4. Mean- field and generalized RPA 

Here we employ a mean-field approximation of the self-energy and a generalized random- 
phase approximation (RPA) of the effective interactions. This represents a useful first 
step toward the fRG treatment, as it allows us to assess the basic behavior and the 
frequency- and momentum-dependence of normal and anomalous vertices. In particular, 
we find that anomalous 3-1-1 vertices have a sizable feedback on certain effective 
interactions. 

4-1- Gap equation 

The system's self-energy can be calculated approximatively by resumming a subset 
of perturbation theory diagrams as in figure O This subset contains all diagrams 




Figure 2. Perturbation expansion for the self-energy. Small circles represent bare 
interactions, lines with arrows represent bare propagators, bold lines with arrows 
represent full propagators, and filled circles represent the self-energy. 

corresponding to bubbles, sunrise diagrams, and trees of either or both. It is 
equivalent to self-consistent Hartree-Fock theory. If we were to restrict the Hubbard 
Hamiltonian to only include the Cooper channel and forward scattering processes with 
zero momentum transfer, the resulting diagrammatic self-consistency equation would be 
exact. Analogous statements can be made about the one- loop fRG treatment for such a 
restricted Hamiltonian [35]. The self-energy equation can be expanded into two distinct 
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parts, 
E 

A- A, 



A(fc) 



/ !2 



E{ky' 

The Matsubara sums are performed analytically. The result is 



Ek 



tanh(/3Efc/2) 



A(r.^ A , rr ^ A(fc) tanh(/3^fc/2) 



(13) 
(14) 



In this approximation, both S and A are constant in momentum space if the external 
field is. The Hartree contribution E can be absorbed into the chemical potential and is 
subsequently ignored. 



4-2. Bethe-Salpeter equation 

A subset of all diagrams of the perturbation expansion for the effective interaction can 
be resummed into a Bethe-Salpeter equation as shown in figure [3l This corresponds 




Figure 3. Bcthc-Salpctcr equation for tlie four-point function. Filled circles represent 
effective four-point functions, empty circles represent bare four-point functions, lines 
which have both ends terminating in a vertex represent full propagators. Note that 
the first two indices refer to incoming lines, while the second pair refers to outgoing 
lines. 



to a generalization of the random-phase approximation in the sense that particle- 
particle ladders in the z/;-representation are included in addition to the usual particle-hole 
bubbles. Symbolically, this Bethe-Salpeter equation reads 

U(1234) = Uo(1234) - ^ U^, MiNs^hiPl, k,P^)^Q MiN2M-iNi{k +^4—^2,^2,^) 

k 

M1...M4, 

X GMiMiik + P4 - P2)GM3M2{k), (15) 

where, as a parameter, i = 1 ... 4 is short for PiNi. By the introduction of the loop 
matrix 



1 , A/4 A/3 



(16) 
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and by exploiting momentum and energy conservation Pi + P2 = Ps + Pa while setting 
p = (p, iUm) := pi — P3, f|T5|) can be rewritten as 



k 

where the Nambu indices of U and Uq are organized into pairs as in (ITTl) . U(pi . . .p^) 
only depends on p in this approximation, as can be seen by solving for U, which yields 



4-3. Numerical setup 

The Matsubara sums in (fTSl) can be performed analytically. The integration over the 
Brillouin zone and the matrix inversion on the right-hand side of (ITSl) is accomplished 
numerically. In the following, a system with t' = t/6, /i = — 1.5t, Uq = 2t, and 
T = 0.035t is assumed. This choice of parameters corresponds to an approximately 
quarter-filled system far below its mean-field critical temperature. A is chosen to be 
real. We mostly employ an external field Acxt = This keeps the numerics well- 

behaved while having no appreciable influence on l]{p) except near p = 0. 

4.4. 3 + 1 effective interactions 

The f2j are generated by diagrams as in flgure [H In the reduced mean-fleld models 
studied previously these vertices remain zero due to the restricted momentum-structure 
of the bare interaction where only particle pairs with total momentum equal to zero 
interact. However, in the general case with scattering of particle pairs with arbitrary 
total momentum they have to be considered. Due to the appearance of an anomalous 
propagator flgure [H these diagrams and hence the are non-zero only in the symmetry- 
broken phase. To obtain an overview of the momentum and frequency dependence of the 
3-1-1 effective interactions, we study the graphs in flgure |H The right-hand graphs show 
that the 3 + 1 effective interactions exhibit no divergence in the Pi — Ps = channel. 
The fij at larger frequencies are suppressed by decreasing the external fleld. From the 
flgure's left-hand diagram, we learn that the f2j's moduli are large around p = and 
around p = (vr, vr) for = 0. 

4.5. Large interaction vertices and Goldstone boson 

The instability toward superconductivity gives rise to a growth of certain combinations 
of the interaction vertices. From previous fRG studies on reduced mean-fleld models 
[26| [27] and from the solution of the Bethe-Salpeter equation, the mean-fleld picture 
of the leading components for zero total incoming momentum and frequency has been 
obtained as follows. There is one strongly growing combination of effective interactions 
that occurs in the flow equation for the anomalous self-energy. It typically reaches 
its maximum at the critical scale or temperature and does not diverge below. 




(17) 
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Figure 4. ^2 for zero energy transfer (left, Aext = 10~^i) and zero momentum 
transfer (right, various Aoxt), approximately quarter filling, Uq ~ 2t, fi = — 1.5i, 
t' — t/Q, T « Tc- r2i,3,4 are similar. 



In the BCS mean-field model, this combination is given hy V + W for real A |26j . 
Another combination of the effective interactions, V — W, diverges also below T^. This 
divergence is associated with the massless phase mode dubbed Goldstone boson present 
if a continuous symmetry, in our case a global ?7(l)-symmetry, is broken. For cases where 
no continuous symmetry is broken as for the half-filled charge-density-wave mean-field 
model with a two-fold degenerate ground state [27] , the cutoff-dependence of the effective 
interaction is comparable to the cutoff dependence of the combination -|- of normal 
and anomalous effective interactions of a BCS mean-field model. 

For the Hubbard interaction studied in this paper, the decay of these large vertices 
away from zero-total- momentum or -frequency becomes important, as other channels not 
present in the mean-field models might pick up large contributions from the symmetry- 
breaking channel. In this section, we investigate V + W and V — W with regard to 
their momentum and energy dependence as well as to the impact of the 3 + 1 effective 
interactions. 



4.5.1. Momentum dependence We analyze the momentum dependence oi V — W and 
V + W at Um = 0. We fit 

^ +7, (19) 



f3\p\^ + A 

to data obtained from numerical calculations for small up to intermediate p, see 
figureO For the non-divergent part V+W, the fit reproduces the numerically-calculated 
behavior. For small momenta, the agreement is better, while for larger momenta, 
deviations become larger, but vanish along the axes. This is due to an anisotropy 
in the numerical data which is not captured by the approximation f[T^ . 

For the divergent part V — W, we note that A as obtained from the fitting 
procedure is proportional to Aext- This is confirmed by calculations for Aext = W~H 
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Figure 5. Left: Blue lines represent fitted to —{V + W) obtained from numerics 
(red crosses). Both plotted in units of A « 0.37t"\ /3 « 1.32t-\ 7 « 3.32t. Right: 
Blue lines represent (fT9)) fitted to — F obtained from numerics (red crosses). Both 
plotted in units of A « l.St-^Acxt, P ~ 0.81i~\ 7 « 3.3M. Both diagrams: C/q = 2t, 
fi = -1.5t, t' = t/6, T « Te, Aext = lO-^i. 



and Acxt = 10~^t. For intermediate momenta along the diagonal (not shown), the fit is 
not optimal as deviations of up to 25% occur. The situation is better for large momenta, 
where deviations of about 10% are common. The peak at p = is well-captured by the 
approximation f|T9l) . which permits us to confirm the dispersion relation for the low-lying 
Goldstone mode in the next section. 

We consider the impact of the 3-1-1 effective interactions. By setting loops combining 
an anomalous and a normal propagator to zero in (ITS!) . the 3-1-1 effective interactions 
are set to zero. The modulus 

Tw/o 3+1 Twith 3+1 rr,n\ 

K := (20) 

J- w/o 3+1 

of the relative change of a given effective interaction combination F upon switching on 
the 3-1-1 effective interactions is plotted for T = V + W and T = V — W in figure [61 We 




note that the magnitude of this change for the combination V + W largely follows the 
magnitude of the 3-1-1 effective interactions in the Brillouin zone. It varies between 
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and 10%. The combination V — W corresponding to the system's phase mode does not 
change to numerical accuracy (note the scale on the z-axis). 

4-5.2. Energy dependence The energy dependence oi V + W and V — is to be 
included in the approximate formula f|T9l) by adding Pyv"^ to the denominator, i.e. by 
considering the approximation function 

^ +7- (21) 



a 



We restrict ourselves to the real part of the effective interactions for brevity. 
Contemplating the left-hand graph of figure [71 we become aware of a shoulder close 




0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 
V [units of t] 



0.04 0.06 0.08 
V [units of t[ 



Figure 7. Matsubara-frcquency dependence oi — {V + W) (left) and W ~V (right) 
calculated far below T^. (green crosses: low resolution, red crosses: high resolution). 
The blue lines are least-squares fits to (|21|) . The parameters used in the fits arc. from 



left to right, {(3y « 0.47t-3,a w 0.31i"\7 w 2.49i}, {a w -0.58t-\ (3^ « 0.46t-3, 



7 ^ 



-3.33t}. Both diagrams: C/q = 2t, ^ ^ -1.5t, t' ^ t/6, T « Tc, Aoxt = 10""* 



to zero frequency, which is not captured by (12T1) . This implies a small non-linearity of 
the dispersion relation of the amplitude mode. 

The ansatz fl2Tl) reproduces the numerically-obtained results for V — W upon 
properly adjusting the parameters, as the right-hand part of figure [7] illustrates. The 
agreement is good for small Matsubara frequencies. In the inset, we see that the large- 
frequency behavior is not captured optimally, which is due to our choice of 7, having 
taken it from the fit in figure [5l However, the exponent 2 for Um is confirmed by the fit. 
This together with the momentum-dependence fit in section 14.5.11 confirms the linearity 
of the dispersion relation of the Goldstone mode connected to — W. 

To understand the influence of the Qi {i = 1 ... 4) on the energy dependence of 
V + W and V — W, we plot the relative change at p = as in ( j20l) . see figure M We 
glean that for V + W, the only appreciable change occurs at = and is quite small 
in the limit Acxt 0. For V — W, the converse is true: Only at = is there 
no change, while even for small z/^ > 0, a sizable change is observed which is only 
suppressed at large Um- This implies that it would not be prudent to neglect the O,, 
while the smallness of their influence on the Um = p = part of -|- TV suggests that 
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Figure 8. Relative change of the energy dependence of V + W (left) and V — W (right) 
upon considering the fli. Both diagrams: Uq = 2t, ^ = — 1.5t, t' = t/6, T << Tc, 



their quantitative impact on the superconducting gap obtained by an fRG calculation is 
likely small. The latter supposition is motivated by the observation that V + W drives 
the flow of the gap, as already shown for the BCS mean-fleld model in [26] . 

5. Functional renormalization group 

In this section, we describe the fRG treatment of the attractive Hubbard model on the 
two-dimensional square lattice at zero temperature. Our main goals are to study the 
flows produced as well as to obtain reliable quantitative results for the order parameter. 
Compared to the partial resummation in the last section, the fRG takes into account 
all one-loop diagrams. The consequential increase in numerical complexity is balanced 
by adopting a cruder discretization of momentum and frequency space. By taking into 
account all one-loop diagrams, the pairing channel is renormalized by other processes 
with non-zero total momentum such as spin- and charge fluctuations. We will flnd 
that this causes a reduction of the superconducting order parameter with respect to the 
mean-fleld result. Furthermore, the fRG treatment allows for a fc-dependence of the 
self-energy around the Fermi surface. 

5.1. Momentum- shell and interaction-flow procedures 
We introduce a cutoff A into our action by substituting 

Go ^ x(A)Go. (22) 

We call 

x(A) = e(|a|-A) (23) 
momentum-shell cutoff function P and 

X(A) = VA (24) 
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interaction- flow cutoff function [12]. Tlie former name is reflected in the expansion of 
tfie system by an infinitesimal shell in momentum space if A ^ A + dA. The latter name 
can be understood by scaling the fields in ^ to \^tp, causing the non- interacting 
part of the action to be cutoff-independent and the interacting, quartic part to be 
proportional to A. For the momentum-shell flow, setting A = Ai > maxfc provides a 
solvable, non- interacting starting point. For A = Af = 0, the original, fully interacting 
system is recovered. For the interaction flow, A = A; = suppresses the interaction 
according to the rescaling of the fields outlined above and provides a solvable starting 
point. A = Af = 1 recovers the original, fully interacting system. 

Employing the momentum-shell cutoff, we introduce a small external field Acxt to 
select the symmetry-broken phase as the endpoint of the fRG calculation. This means 
our results in this case apply to a situation with spontaneously-broken symmetry only in 
the limit Aext — ^ 0. For the interaction-flow cutoff, we add zero in the form Acxt — Acxt 
to the quadratic part of the action, including Acxt in the non-interacting part and —Acxt 
in the interacting part of the action. This approach has been shown to yield symmetry- 
broken results at zero external field [28] for A = Af. 

5.2. fRG equations in the Katanin truncation 

The flow equations which we employ are shown in figure [9l where the single-scale 




Figure 9. IPI flow equations in the Katanin truncation. Internal lines stand for full 
propagators, a slashed full line represents the single-scale propagator §, a disc with two 
legs stands for the self-energy, a disc with four legs represents the effective interaction, 
a dot denotes differentiation with respect to the cutoff A. 

propagator S = — G(?a (Gq^) G is used. The name can be understood by noting that 
because of the differentiation acting on 0, § lives only on the energy scale A for the 
momentum-shell flow. These flow equations can be obtained by truncating according 
to Katanin [25] the infinite hierarchy of one-particle irreducible (IPI) flow equations 
of [11]. Solving them yields the exact order parameter and effective interaction in the 
thermodjTiamic limit for any model with only two-particle interactions which is mean- 
field exact [23 [361 [26]. 
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Symbolically, the flow equation for U reads 

11(1234) = lJNiN2]VhM2iPl,P2, k)UM3MiN3Niik, -k + p, P3)GM3Mi{k)GMiM2i-k + p) 

k,M 

+ ^ UMi7V2M4Af4(^ - 9,P2, k)l]NiM3N3M2iPl^ k,p-i)dp, M^Auik)'^ NhAhik - q)] 
k,M 

- "^M^N^MiNiik - q',Pl, k)l]N2M3N3M2{P2, k,ps)dA [GM3M4 (^)Ga/i A/2 " q')] > (25) 
k,M 

where p = pi +P2, Q = Ps —Pi, and q' = ps —p2- Note that the last two terms correspond 
to the second diagram pictured above, as the antisymmetrization is denoted explicitly 
in the symbolic equation. Prefactors arising from Fourier transformations are absorbed 
into the summation symbols. 

The symbolic form of the flow equation for the self-energy reads 

^NiN2{p) = ^MiN2M2NAk,P,k)x{k)^^^^{k). (26) 

k,MiM2 ^ 

Here, the single-scale propagator is rewritten into a form better suited for numerical 
evaluation [27] . 

5.3. Numerical results 

In the following, we simplify (125|) and (l26l) and subsequently present results stemming 
from a numerical quadrature of the simplified flow equations. While we focus on the 
momentum and frequency structure of the effective interaction in section HI in this 
section we are mainly interested in determining the order parameter beyond mean-field 
theory, taking into account all one-loop diagrams in the fRG. However, this comes at 
the price of an increase in numerical complexity, which we balance by adopting a cruder 
discretization of the effective interaction. We first drop the frequency-dependence of 
U. Consequently, no energy dependence of S is generated in the flow according to 
f l26|) . We can therefore do all Matsubara sums analytically. We furthermore drop the 
normal part of the self-energy, thus neglecting Fermi surface shifts. All quantities can be 
considered real in our approximation. A parametrization taking into account partially 
the momentum dependence of the effective interaction is the so-called "patching" [91137]. 
The Brillouin zone is split into patches as shown in figure [TOl The functional value of the 




Figure 10. Patching of the BriUouin zone used with varying total number of patches 
in the numerical calculations. Every patch covers the same angle. 
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effective interaction at the mid-angle of the patch and on the Fermi surface is assumed 
to be also valid for all other momenta in the patch. Thereby, only the momentum 
dependence along the Fermi surface is taken into account. This is motivated by the 
finding that the dependence perpendicular to the Fermi surface is irrelevant in the RG 
sense for the momentum-shell flow without symmetry-breaking ([5], [6]). We specify 
the momentum arguments of U as triples of patch indices, so \^(1,2,3) refers to any 
^(Pi;P2;P3) with pi in patch 1, p2 in patch 2, and pa in patch 3. Note that to obtain 
momenta in the ip representation from these patch numbers, signs must be adjusted 
according to (Q. 

5.3.1. Momentum- shell flows We replace the step function from fl23l) by 

^■('\) = 1 - ?=TT- 

This simplifies the numerical calculation significantly, although one could suspect that 
the sharp cutoff might eliminate an integration. However, in the Katanin truncation, 
terms proportional to the derivative of the order parameter arise in addition to the 
terms proportional to the derivative of the cutoff function. Note that our results do not 
depend on a in a reasonable range, which adjusts the sharpness of the cutoff. 

In figure [TTl we plot flows of the order parameter in the first quadrant of the 

0.06 

0.05 

^ 0.04 
"o 

.3 0.03 
< 

0.02 
0.01 


le-05 le-04 0.001 0.01 0.1 1 10 le-05 le-04 0.001 0.01 0.1 1 10 
A [units of t] A [units of t] 

Figure 11. Momcntum-shcU flow of the order parameter. Twelve patches, quarter 
mihig: ^ = -1.41t, t' = -O.li; C/q = 1.5t, Aoxt lO^^t (left), Acxt = 5 x 10"''* 
(right). 

Brillouin zone. We note that the flow saturates for the two external field strengths 
employed. However, if we employ too small an external field, the effective interaction 
flows diverge and the method breaks down. Fortunately, the external field can be chosen 
at least a hundred times smaller than the final value of the superconducting gap. In 
calculations for mean-field models [27], the modification of the order parameter due to 
the external field is less than 20% in this case. 

From the resummation studies in section HI we have learned that the Goldstone 
phase mode V — W dominates for zero momentum transfer in the Nambu formalism, 
which is equivalent to zero total momentum if we transform back to the original ip fields 
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according to (jH])- In figure W2\ we plot several fiows from this channel. They dominate in 




Figure 12. Momentum-shell flows of the effective interaction V — W associated 
with the Goldstone mode. Twelve patches, quarter filling: fj, = — 1.41t, t' = — O.lt; 
Uq = 1.5t, Aext ~ IQ^^t (left), Aext = 5 X 10~^t (right), numbers in parentheses are 
patch numbers to be understood within the Nambu formalism. 



magnitude all other effective interaction flows, including those not shown in the flgure, 
and exhibit a spread of ~ 10% under variation of the momentum combination. They 
are quite similar to the flows found in [26]. Note that the maximal flnal value modulus 
increases stronger than proportionally to A^^- This is a precursor of the divergence to 
arise for even smaller Aext- 

Complementing the phase mode is the effective interaction combination V + W 
driving the order parameter. Flows of + are plotted in flgure [131 The spread of the 




le-05 le-04 0.001 0.01 0.1 1 10 le-05 le-04 0.001 0.01 0.1 1 10 



Figure 13. Momentum-shell flows of the effective interaction V + W associated 
with the amplitude mode. Twelve patches, quarter filling: /i = —lAlt, t' = —O.lt; 
Uo = 1.5t, Aext = I0~^t (left), Aext = 5 x 10~'^t (right), numbers in parentheses are 
patch numbers to be understood within the Nambu formalism. 



flnal values increases with decreasing Aoxt, while the average remains approximately 
constant. The graphs resemble those found when a discrete symmetry is broken 
|27] . However, the monotonicity below the critical scale which had been established 
analytically for the mean-fleld model is violated by one of the flows in the right-hand 



Superconductivity in the attractive Hubbard model: fRG 



17 



figure. Overall, all momentum- shell flows exhibited resemble the flows for the mean-field 
models. 

5.3.2. Interaction flows In this section, we employ the interaction flow cutoff ( !24l) 
with an initial anomalous self-energy As,e.,i = Aext and a counterterm Ac = Acxt which 
we include in the bare propagator. In figure [Ml we show examples of the flow of the 




I ' ' ' ' 1 0.031 I ' ' ' ' 1 

0.2 0.4 0.6 0.8 1 0.95 0.96 0.97 0.98 0.99 1 

A A 



Figure 14. Interaction flows of A = Ac — x^s.e. for patcii 1 at C/p = 1.5t, t' = —O.lt, 
jjL ~ —I.Alt (quarter filling), Ac ranges from 0.05i to 0.19^ and can be read off at the 
intersection of flow line and A-axis. 

superconducting order parameter A = Ac — x^s.c. for various counterterms. We note 
that if we could treat the flow equations without truncation and further approximations, 
all counterterms should produce the same final self-energy and vertices. Indeed, the 
counterterm-fRG for a half-filled charge- density- wave mean-field model exhibited a 
single strong attractor for the flow of the self-energy 

For the more general model studied here, the truncation error becomes noticeable, 
and for certain values of the pairing counterterm, the interaction vertices still diverge. 
The flows are terminated if the maximal effective interaction exceeds four times the 
bandwidth, i.e., if the flow leaves the weak-coupling regime, or if A reaches Af = 1. We 
see in the left diagram that the flows no longer converge onto a single strong attractor 
as for the mean-field case [28]. Studying the final 5% of the flow shown in the right 
plot of figure [TH reveals that there exists a set of small counterterms for which the 
flows cross and leave the weak-coupling regime, and a minimal counterterm for which 
the effective interactions remain within the above limit up to Af = 1. The flows for 
counterterms whose strength is just beyond this minimal counterterm terminate closer 
to each other than do the flows for larger counterterms. The minimal counterterm is 
approximately twice as strong as the final order parameter value obtained in the minimal 
counterterm flow. The behavior described above is interpreted as being analogous to the 
strong-attractor behavior in the mean-field case. The overestimation of the final order 
parameter value for larger counterterms is attributed to the overestimation of the order 
parameter in the flow due to the neglect of the order parameter's energy dependence and 
momentum dependence perpendicular to the Fermi surface. Consequently, we choose the 
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minimal order parameter obtainable by terminating the flow at Af in the weak-coupling 
regime as our approximation for the physical order parameter. Another way to argue 
for the so chosen counterterm is that the continuous symmetry breaking necessitates a 
divergent vertex exactly at the end of the flow in order to render the Goldstone mode 
correctly. This inflnite value is suppressed by our approximation, but the flow for 
the chosen counterterm yields the maximal flnal value for the Goldstone mode within 
the above limit. Finally, note that the order parameter obtained by these rules is 
approximately 20% smaller than the value obtained with the momentum- shell method 
in section 15.3.11 which is enhanced by a flnite, albeit small, external fleld. 

5. 3. 3. Order parameters flgure [15] shows the magnitude of the superconducting order 



'S 
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-2.5 -2 -1.5 -1 -0.5 0.5 



)i [units of t] 

Figure 15. Magnitude of the superconducting order parameter A as calculated via the 
interaction flow procedure for varying interaction strength Uq and chemical potential 
H, t' = -0.lt 

parameter A versus the chemical potential as calculated using the interaction flow 
method including a counterterm as described in the previous section. For small bare 
coupling f/o, A is maximal if the van Hove points lie on the bare Fermi surface, which is 
the case for = 4t'. For larger bare coupling, the maximum gradually shifts to smaller 
/i. We compare the results from the fRG to results from mean-fleld theory in flgure [T6l 
We note that the suppression of the mean-fleld values by the fluctuations taken into 
account by the fRG varies between a factor 1.6 and 2.0. This is quite similar to values 
found by second-order perturbation theory [33]. Table [1] compares our results to results 
from a combination of fRG and mean-fleld theory [38] for flnite t' . While the values are 
of the same order of magnitude, the results obtained by employing only the fRG are 
larger by between 20 — 30%. 
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Figure 16. Dependence of the superconducting gap on tiie effective interaction as 
calculated via the interaction flow procedure, quarter filling: t' = —O.lt, jj, = —lAlt 
(left), half filling:t' = 0, ^ = (right). For the quarter-filled case, fits to a exp(— /3/f/o) 
are included, where a and f3 are real fit parameters. 



Table 1. Superconducting order parameter A in units of t for t' = —O.lt, 
/i = — 0.2612t (half filling), comparing fRG results to results from a combination of 
fRG and MF. *: value obtained using linear extrapolation. 



U[t] 


fRG 


fRG and MF 


1.0 


0.02 


0.02 


1.5 


0.11 


0.07 


2.0 


0.22 


0.15 


2.5 


0.35 


0.26 


3.0 


0.49 


0.38* 



6. Conclusions 

We have applied the fermionic one-particle irreducible (IPI) version of the functional 
renormalization group (fRG) in the Katanin truncation to the symmetry-broken 
phase of a two-dimensional, non-mean-field model exhibiting superconductivity at zero 
temperature. The attractive Hubbard model has been chosen as the prototype system 
because it is sufficiently well-studied to provide a good testbed for new many-body 
methods. In this case, particle-hole as well as particle-particle loops are important for 
the symmetry-breaking properties of the system. The fRG takes both into account on 
an equal footing. 

Our analysis has proceeded in three steps. First, the usual Nambu formalism has 
been extended to cover 3-1-1 anomalous effective interactions featuring an odd number 
of incoming legs, whose behavior and importance was a focus of our interest. We have 
found that for the case of superconductivity, the effective interaction had to be extended 
into a 4 X 4 matrix containing certain redundancies due to symmetries. This extension 
followed from a substitution of the fields in which the model is originally formulated by 
Nambu fields which are chosen specifically to produce an action of canonical shape while 
encompassing singlet Cooper pairs in the quadratic part. 3 + 1 effective interactions were 
found to arise naturally within this formalism. 
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In the second step, a study of the system has been performed by means of a 
resummation of a subset of all perturbation theory diagrams. The strength of the 
symmetry breaking was determined by a mean-field gap equation. A class of diagrams 
containing bubble chains, ladders, and all chains comprising both was resummed into 
a Bethe-Salpeter equation. In the context of this equation, the impact of the 3 + 1 
effective interactions was studied by comparing effective interaction values calculated 
with and without the 3 + 1 effective interactions. This impact was found to be zero for 
the zero-frequency part of the phase mode, and to be dependent on the external field 
for non-zero frequency for both phase and amplitude mode. 

In the third step, the fRG has been numerically applied to the attractive Hubbard 
model. Two distinct variants of the fRG were employed. First, a momentum-shell 
cutoff was used together with an external field providing a finite initial order parameter. 
The Goldstone mode divergence is regularized by the field. Despite the approximate 
nature of the flow and the explosive character of the fRG differential equation, flows 
reaching down to zero scale were found employing external field strengths two orders 
of magnitude smaller than the final order parameter values. The relative error of the 
results for the order parameter due to such an external field was estimated to be at 
most 20%, judging from earlier calculations for mean-field models. For even smaller 
external field values, a regime of diverging flows was found. Second, we have calculated 
interaction flows, balancing an external field with a counterterm increasing continuously 
with the flow. We have noted that vestiges of strong-attractor behavior found previously 
for a mean-field model are also present for the attractive Hubbard model flows. They 
permit the determination of order parameter values which are comparable to results 
from the literature. While these values proved to be smaller than values calculated 
by mean-field theory alone, as expected, they are generally slightly larger than results 
obtained by combining symmetric-phase Wick-ordered fRG calculations and mean-field 
theory. This deviation is likely due to small differences between the IPI and Wick- 
ordered fRG schemes, since symmetric-phase IPI calculations already yield larger critical 
scales than the corresponding Wick-ordered calculations. In summary, the IPI fRG in 
the Katanin truncation was found to yield fairly accurate quantitative results for weak 
and moderate interaction strengths (up to Uo/t = 3). We reckon that the diverging 
flows present for very small external fields in our calculations can be pushed back by 
improving the resolution of the discretization in both frequency and momentum space. 
This would furthermore improve the accuracy of the results. The method may also 
be employed to study symmetry breaking in interacting-fermion models for which a 
stronger momentum dependence of the order parameter is expected, for example to 
study the repulsive Hubbard model. The method can also be applied to situations 
with competing order. The interaction-flow scheme with possible symmetry-breaking 
in two different channels has recently been tested in a one-dimensional model with 
competing long-range correlations. There it allowed for a correct determination of a 
quantum critical point between a charge-gapped phase and a compressible phase with 
dominant superconducting correlations [39] . Finally, it will also be interesting to analyze 
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how non-zero temperature affects the renormahzation of the mean-field-hke instabihty 
encountered in the ffows. 
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